function y = prb_inv (fun, x, p, q, x0)

   ## usage:  y = prb_inv (fun, x, p, q, x0)
   ##
   ## inv. probit transform using fun to fit cdf

   y = x ;
   I = isnan(x) ;
   w = normcdf(x(!I)) ;
   w = (w-q) ./ (1-q) ;
   p = mat2cell(p(:)', 1, ones(1, length(p))) ;
   w = feval([fun "inv"], w, p{2:end}) ;
   II = !isnan(w) ;
   w(II) = max(x0, w(II)*p{1} + x0) ;
#   w(!isfinite(w)) = inf ;
   w(!II) = 0 ;
   y(!I) = w ; y(I) = nan ;

endfunction
